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We examine the effects of turbulence on elliptic instability of rotating stratified 
incompressible flows, in the context of the Lagrangian-averaged Euler-Boussinesq- 
alpha, or LAEB— a, model of turbulence. We find that the LAEB— a model alters the 
instability in a variety of ways for fixed Rossby number and Brunt- Vaisala frequency. 
First, it alters the location of the instability domains in the (7, cos 0)— parameter 
plane, where 9 is the angle of incidence the Kelvin wave makes with the axis of 
rotation and 7 is the eccentricity of the elliptic flow, as well as the size of the asso- 
ciated Lyapunov exponent. Second, the model shrinks the width of one instability 
band while simultaneously increasing another. Third, the model introduces bands 
of unstable eccentric flows when the Kelvin wave is two-dimensional. We introduce 
two similarity variables-one is a ratio of the Brunt- Vaisala frequency to the model 
parameter Tq = 1 + a^/?^, and the other is the ratio of the adjusted inverse Rossby 
number to the same model parameter. Here, a is the turbulence correlation length, 
and P is the Kelvin wave number. We show that by adjusting the Rossby number and 
Brunt- Vaisala frequency so that the similarity variables remain constant for a given 
value of Tq, turbulence has little effect on elliptic instability for small eccentricities 
(I7I ^ 1). For moderate and large eccentricities, however, we see drastic changes of 
the unstable Arnold tongues due to the LAEB— a model. Additionally, we find that 
introducing anisotropy in the vertical component of the transported velocity field 
merely alters the definition of the model parameter Tq, which effectively reduces the 
original parameter value. When the similarity variables are viewed with the new 
definition, the results are similar to those for the isotropic case. 
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PACS numbers: 47.20. Cq, 47.20.Ft, 47.27.Eq, 47.27.Cc 
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I. INTRODUCTION 

This paper is the third in a continuing investigation of the effects of a turbulence clo- 
sure model on the classic solution of elliptic instability. The firsiji examined the effects of 
the Lagrangian-averaged Navier-Stokes-alpha, or LANS— a, turbulence closure model in an 
inertial frame. The second^ examined the effects of this turbulence model in a rotating 
frame, as well as the effects of a few related turbulence closure models. In the present paper, 
we continue the investigation by examining the effect of stratification with rotation in the 
inviscid LANS— a model, also called the Lagrangian-averaged Euler-Boussinesq-alpha, or 
LAEB-a, model.-^ 

Turbulence models are ubiquitous in fluid mechanics today. They are used with the hopes 
that they adequately describe the evolution of real-world flows, governed by the Navier- 
Stokes equations. The question that arises, then, is how do individual turbulence models 
affect classical solutions of the Navier-Stokes equations? A natural thought process is that 
a model which preserves classical solutions is preferred over a model which drastically alters 
the solutions. However, imposing such a strict requirement on a turbulence model is unre- 
alistic. An alternative thought process asks how to change the problem so that the effects 
of the model on classical solutions is negated while using the low-resolution advantage the 
turbulence model possesses. Our present investigation has two goals. First, we seek to un- 
derstand how the LAEB— a turbulence model affects the classical solution known as elliptic 
instability. The model is formulated in a way that preserves the existence of the solution, 
something which is not possible with every turbulence model. We find that the model alters 
the instability, which leads to our second goal. Namely, we seek a reformulation of the 
problem so that the effects are negated. Turbulence models are being used in applications 
such as ocean circulation. Our motivation is to determine how turbulence closure models 
affect instabilities before implementing the models in such applications. 

The LAEB— a model is a turbulence closure model which has two velocities, a transport 
velocity and a transported velocity. The transport velocity u is smoother that the velocity 
V (or momentum) which it transports, given by the relation v = (1 — a^A)u in the isotropic 
version. (An anisotropic version will be given in Section 5.) This approach was first intro- 
duced by Leray^ as a regularization of the Navier-Stokes equations. Here, the parameter a 
is interpreted as the turbulence correlation length. The idea is that in numerical simulation. 
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a is set to be the grid size, and that the dynamics which occur at scales smaller than a are 
swept along to scales larger than a by the transport velocity u. 

Elliptic instability is a mechanism by which two-dimensional swirling flows become un- 
stable to three-dimensional perturbations. The instability in the Euler equations was in- 
vestigated by Lord Kelvin^ for the case of circular streamlines and by Bayly^i for elliptic 
streamlines. Further work by other o^i^^i^-^i^^i^^i-^'^i'^^i'^^ examined the effects of rotation, strat- 
ification, and magnetic fields on the instability. We refer the reader to a recent review^ for 
details about the history of elliptic instability. It is a member of the Craik-Criminale, or 
CC, class of exact solutions to the Navier-Stokes equations.— The velocity field in this class 
of solutions is the sum of a flow linear in the spatial coordinates and a Kelvin wave. The 
present authors have shown that this class of solutions are valid in a family of turbulence 
closure models.— This investigation is continued here. 

The paper is organized as follows. In Section 2, we state the LAEB— a equations and 
discuss the CC class of solutions for these equations. We examine a specific member of this 
class in Section 3, namely elliptic instability. We combine analysis and numerical simulation 
to give a full description of elliptic instability in the isotropic LAEB— a model in Section 4. 
The analysis yields two similarity variables which fully describe the behavior of the instability 
for elliptic flows of small eccentricity. We briefly discuss the effects of the anisotropic version 
of the LAEB— a model in Section 5. We show that the effects of anisotropy are similar to the 
isotropic version, once the similarity variables of Section 4 are redefined. Finally, Section 6 
concludes the paper with a summary of our findings, including a conjecture as to why the 
model alters the unstable Arnold tongues for moderate and large eccentricities. 

II. CC SOLUTIONS IN LAEB-a EQUATIONS 

The isotropic version of the LAEB— a equations^*^ are 

(9tv + u ■ Vv + (Vu)^ ■ V + 2r2 X u + hgQ^ + V (^p - -|up - — | Vu 

+ u ■ V6 = 0, 
divu = 0. 

The variables are as follows: 

• V = (1 — q;^A)u is the fluid velocity field. 



= F, (1) 

(2) 
(3) 
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• u is the mean transport fluid velocity field, 

• p is the pressure, 

• b = b{z) is the dimensionless buoyancy term which reflects the deviation from mean 
density [that is, we have taken p = po(l + b)], 

• Q = is the Coriolis force, 

• g is the gravitational constant, 

• F is the contribution of all external body forces, 

• in index notation, (Vu)ij = djUi = Uij. 

The equations reduce to the classic Euler-Boussinesq (EB) equations upon setting a = 0. 
The equations admit solutions (u,p, 6) = (Uo,-Po)^o) given by Uq = S(t) ■ x + U(t), Pq = 
■ T{t) ■ X + P(t), and bo{z) = s{h — z), provided 

dtSij + SirSrj + 2^lei3pSpj - sg63i63j = Mij, (4) 

S^^ = 0, (5) 

where we sum over repeated indices. Here, Vq = (1 — a^A)Uo = Uq, and the (symmetric) 
Hessian Ai is defined by = —didjF, where 

F = p- j F ■ rfx+ {dtUg - SgrUr + shg63g)xg . 

The Brunt- Vaisala frequency, defined for the EB equations as A^^ = —{g/po){dp/dz), is 
= sg with s > 0. 

One can construct a new exact solution of the LAEB— a equations of the form 

u = Uo + Ui, b = bo + bi, p = Po + Pi. (6) 

Consequently, v = Vq + Vi, where Vi = (1 — q;^A)Ui. The generic equations of motion for 
Vi, bi, and Pi are 
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dtYi + Uo ■ Wi + Ui ■ Wo + Ui ■ Wi + (VUo)^ ■ Vi + (VUi)^ ■ Vo 

+ (VUi)^ ■ Vi + 2^2 X Ui + 6i(?e, - (VUq)^ ■ Ui - (VUi)^ ■ Uq (7) 

+ V (Pi - i|Uip - y iVUip - a2(vUo).,(VUi),, j = 0, 

dM + Uo ■ V61 + Ui ■ V60 + Ui ■ V61 = 0, (8) 
divUi=0. (9) 

Equations ((Tj) and (jH)) contain two nonlinear terms. We choose the disturbance for the 
buoyancy and transport velocity field u in the form of a traveling Kelvin wave: 

Ui = /ia(t)e^^'^, hi = sfib{t)e"^^, Pi = fipi{t)e'f^^' + fi%{t)e^''^'^ , (10) 

where ip{'x,t) = k"^(t) ■ x + S{t), (3 is the wave number chosen so that |k(0)| = 1, and /i is 
a scaling parameter so that |a(0)| = 1. Consequently, the transported velocity field has the 
form Vi = /iT(i)a(t)e*^'^, where 

T{t) = l + a'(3^\\i{t)\^. (11) 

We refer to solutions of the form © with (jlOp as Craik-Criminale, or CC, solutions.— A CC 
solution can be viewed as the sum of a 'base fiow' Uq plus a 'disturbance' Ui in the form 
of a traveling Kelvin wave. The incompressibility condition (jUj) yields transversality of the 
wave 

k^ ■ a = 0. (12) 

Upon inserting (jlUp into ((Tj) and (jHI), the terms in ((7j) and (jHI) nonlinear in the disturbance 
identically vanish due to wave transversality, implied by (jl2|) . The remaining terms are 
linear and constant in x. The CC solution ansatz requires that the coefficients of these 
terms vanish separately, thereby yielding the following system of equations for the wave 
vector k, wave amplitude a, and pressure harmonics pi,2'- 

dtk^ ■ X + k^ ■ 5 ■ X = 0, (13) 
rft(Ta) + i(3Ta{dt6 + ■tJ)+S^ ■ (Ta) + n x a - {i^pi - a^(3^k^ ■ 5 ■ a)k 

(14) 

+ N^be, = 0, 

sdti) + is(3b{dt6 + k^ ■ U) + (Vbo) ■ a = 0, (15) 
P2 + (T-l)|a|2 = 0. (16) 
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Here, 11 = 212 + curl Uq is the total vorticity of the rotating coordinate system and the base 
flow. Without loss of generality, we set dt5 + k-^ ■ U = 0. We take the gradient of (fT^ and 
obtain the following system of equations: 

rftk + 5^-k = 0, (17) 
dtija) + 5^ ■ (Ta) + H x a - Pk + N%e, = 0, (18) 
dtb-SL-e, = 0, (19) 

subject to p2|l . The term P represents the coefficient of k in |T^ . At this point, one may 
assume that a(t), b{t), ipi{t), and P2{t) are real- valued functions. 

The operator {dt + S'^-) in ()17|) and (fTHj) is the total derivative of a vector field in a 
Galilean frame moving with Uq. Thus, these equations state that the wave vector k is 
frozen in the fiuid, and the scaled amplitude Ta in the Galilean frame moving with Uq is 
affected by the total vorticity of the system, by pressure, by the buoyancy, and by the base 
flow. Once the CC solutions are established, introduction of the factors of T in equations 
fllfi|l and ()18|1 reveals the key differences between the elliptic instability analyses for the EB 
equations and the mean turbulence equations of the LAEB— a model. 

As is common when working with incompressible fluids, the pressure term can be removed 
from the problem by taking the dot product of ()18|1 with k and using the fact that k"^ ■ a is 
an integral of motion: 

k^ ■ (T(5 + 5^) ■ a + n X a + N%e,) 



P 



|2 



|k| 

Once ()17|) is solved for k, the equations for a and b can be written compactly as 





Ar(t;...) I . (20) 



Here, A/" is a 4 x 4 matrix whose coefficients contain elements from S and k, and the dots 
indicate parameter dependence. 

Multi-harmonic disturbances. Finally, we point out that for the EB equations (corre- 
sponding to a = 0), the perturbation can be over a sum of harmonics of e^^^. In this case, dif- 
ferent harmonics do not interact. This phenomenon changes slightly for the LAEB— a equa- 
tions. For example, if we were to consider Ui = jjsS^^ exp^iPip) + jj^a^'^^ exTp{2iPip) and 
similarly for Vi and bi, then the correct form for the pressure perturbation Pi would carry 
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the first four harmonics of exp{if3il)). The pressure coefficients Pi{t) of the first, third and 
four harmonics would be purely real. However, the coefficient of the second harmonic would 
be complex-valued such that the real part balances a contribution from a^^^^ and the imagi- 
nary part from a*^^). Thus, different harmonics do interact, but the interaction is absorbed 
by the pressure. The remainder of the paper considers only single-harmonic disturbances. 



III. ELLIPTIC INSTABILITY 



The classic problem of elliptic instability resides among the CC solutions. The base flow 
for elliptic instability is 



S 



UJ 





1+7 




-1+7 





U = 0. 



(21) 



This base flow is a rigidly rotating column of fluid whose streamlines are ellipses with 
eccentriticy 7 and vorticity 20)6^. We rescale time by uj so that the Coriolis parameter 
r2 is replaced by fi' = fi/o), which we interpret as an inverse Rossby number.— We now 
suppress the prime notation. The extreme values of the solution in (j2H) represent pure shear 
at 7 = ±1 (along different axes) and rigid body rotation at 7 = 0. The solution of ()17|1 is 



^ sin Q cos i ^ 
K sin Q sin t 

cos B 



(22) 



where i = ty\ — 7^, /t = (1 — 7)/(l + 7), and Q is the polar angle the wave vector k 
makes with the axis of rotation e^. We note that a more arbitrary initial orientation of k is 
equivilaent to a shift in the starting time t ^ t — to. 

In general, the solution of ()20|) must be simulated numerically. Since the coefficient matrix 



A/" in fl20j) is periodic with period r = l-nj-^X — 7^, it follows that (PUI) can be rewritten 
as a pair of Schrodinger equations with periodic potentials, shown in the appendix. Thus, 
we analyze the system using Floquet theory One needs to determine the eigenvalues pi of 
the monodromy matrix V{t), where dtV = ■ V with the initial condition 'P(O) = X4 = 
diagjl, 1, 1, 1}. The Lyapunov growth rates, which exist only if maxj \pi\ > 1, are given by 
a = ln(maxi |3fJ(pi) |)/r. 
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IV. ISOTROPIC LAEB-a MODEL 

We begin by examining the isotropic model presented in Section II. The classical re- 
sults for the EB equations were obtained by Miyazaki and Fukumoto^ and Miyazaki.i^ 
Rather than review the classical results separately, we will study the solution in the isotropic 
LAEB— a model directly, regaining the classical results upon setting a = 0, which is equiv- 
alent to setting T(t) = 1. To illustrate the effects of the model, we plot the classical results 
next to those generated by the model. The piece de resistance of the present paper is the 
introduction of the similarity variables 

2|fi + l| 



J-o 



Tn 



(23) 



where 



To = 1 + a^(3\ (24) 
Since a has dimensions of length, it follows that Tq is a dimensionless parameter. 

A. Circular streamlines— exact solutions 

The circular case 7 = can be analyzed analytically. Here, |k(t)| = 1, and any solution 
to ()20|) can be written as 




where 



ai = sin(u;t + 0)k^i + 
a2 = — cos(u;t + </))k_|_i 



as 



1 




a2 \ / ag \ / a4 



(25) 



1 + q^x 
1 



^ q^X cot 6 cos t + qxt sin 
q^cotesint-qxtcost 
1 



cos(ci;t + </))k_|_2, bi 
sin(u;t + 0)k_L2, h 



) 



cos (tut + (f>) 



sin(u;t + (j)) 



a4 = -gxk 



±2, 
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Here, k 



■±i — 



[cos 9 cos t, cos 9 sin t, — sin 9] 




■±2 — 



[— sint, cost, 0]'^, 



2Ccos9y^TT7x 



and q = tan9/(. For a = 0, these solutions are the exact version of the 0{N'^) solutions of 
Kerswell.ii Additionally, for A^^ ^ Q (x = 0), we regain the LANS— a solutions^ for a > 
and the classic Euler solutionsi^ii^ for a = 0. 

To compute the monodromy matrix V{27r), we use the solution in (j25|) and determine the 
coefficients q for each column of the monodromy matrix by the condition V{0) = l^^. The 
resulting eigenvalues are pi^2 = 1,^3,4 = exp(±2ico'7r). Since all of the eigenvalues satisfy 
\pi\ = 1, the circular case (7 = 0) is stable. Furthermore, the angle of critical stability, 
that is, the parameter value at which instabilities will set in as 7 increases from zero, is 
determined by the condition pi = 1, or equivalently, 2ct;7r = ±2n7r, where n = 1,2,3,.... 
This corresponds to 



The first conclusion is that stratification and rotation increases the angle of critical stability 
(decreases cos^). The different values of n indicate the primary, secondary, tertiary, etc., 
instability domains. These different instability domains, or fingers, have been extensively 
studied.— i^^i^^ The principal finger, n = 1, is the main instability domain for elliptic insta- 
bility, and it is call the 'subharmonic frequency.' The secondary finger, n = 2, corresponds 
to a resonance phenomenon with internal gravity waves, and thus is called the 'fundemental 
frequency.' The remaining fingers, n = 3,4, . . ., naturally are referred to as 'superharmonic 
frequencies.' In the absence of statification, only the subharmonic frequency (n = 1) is 
physically interesting from the viewpoint of elliptic instability, although rotation increases 
the width of the other frequencies. (These fingers are displayed in Figs. below.) Since 
the expression for cos^ is independent of To explicitly, the second conclusion we draw is 
that the LAEB— a model has no effect on the circular case, provided that the values of A^^ 
and Q change as a changes so that x ^'^^ C remain fixed. 




(26) 



B. Slightly elliptic streamlines perturbation analyses 



The parameter values which yield critical stability for the circular case will produce 
instability in the elliptic case for small eccentricities. We derive an analytical expression for 
the leading order growth rate of this instability to 0(7^), defined for the LAEB— a model 
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as 

^ ^U^M' + N%'). (27) 

This is a natural extension of the definition given by Kerswel^^ for To = 1, and differs 
shghtly from that in previous work on the LANS— a modeL^^ By taking the dot product of 
(fT^ with a, multiplying by and adding, we find that the following expression for 
cr: 

-a^ ■ ■ (Tpa) 
a = — = -7aia2 

z 1 

valid for all 7, where ai^2 are the first and second components of a, respectively. We insert 
the circular solutions (that is, for 7 = 0) into the above equation, say ai and bi, and average 
over a period of the solutions. Typically,— the average will vanish except when u ± 1 = 0, 
a condition which again yields (j2EI)- A maximum is attained when = ±7r/4, yielding 

We see that for base fiows with small deviations from circular, that is, for fiows with I7I -C 1, 
the growth rates are linear in the eccentricity. Of particular interest is the fact that Cmax is 
a function of all three parameter x, (, and Tq, rather than just the first two. We conclude 
that although properly adjusting the Brunt- Vaisala frequency and Rossby number preserves 
the effects the model has on circular flows, the growth rate in the LAEB— a model is a 
decreasing function of the model parameter Tq. Numerical simulations below will show 
that for moderate and large eccentricities, the effect is more dramatic. In particular, the 
subharmonic flnger (n = 1) has a maximum growth rate of a^^x = 7/(2To), which occurs 
at C = 1 (r2 = — 1/2). For this parameter value, the critical angle is cos^ = 1, that is, when 
the wave vector is parallel to the ellipse's axis of rotation. This is referred as the zero tilting 
vorticityM^ We will expand on this case in a separate section. 

The expression for the critical angle separates the parameter space into four distinct 
regions. We describe each region by its characteristic behavior. We preserve the ( and 
X notation since the behavior will remain unchanged under the redeflnition of these two 
variables for the LAEB— a model in the next section. 

(I) Both the numerator and denominator of the radicand have the same sign and the 
ratio is less than unity. For parameter values in this region, we have cos 6* < 1. 
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Thus, slight perturbations in the columnar vortex's eccentricity from circular will 
yield exponentially growing amplitudes. This is due to the fact that the fingers touch 
the line 7 = at a point (a consequence of parametric resonance and the fact that 
7 = is stable). For a given finger with index n, we have that either cos^ 0"^ as 
X — ^ ('^^)~ and it vanishes for larger values of x or cos^ — > 1~ as x oo? iii which 
case it never vanishes. 

(II) Both the numerator and denominator of the radicand are positive and the ratio is 
greater than unity. Parameter values in this region correspond to cos 6' > 1 and 
experience a window of nonzero eccentricities for which the Kelvin waves are bounded 
(i.e. stability). If we could extend the parameter plane to include values of cos 6 larger 
than unity, then the fingers which fall into region (II) would also experience instability. 
The unique property here is that cos^ +00 as x — * (C^)~- 

(III) Both the numerator and denominator of the radicand are negative and the ratio is 
greater than unity. For these parameter values, the fingers whose index corresponds 
to n < ( have vanished and will not reappear. They enter this region at x = (C^)"*" 
with cos 6' = +00 and has the feature that cos 6' ^ 1^ as x ^ 00. 

(IV) The numerator and denominator have different signs. Parameter values which fall into 
this region are particularly interesting. For these values, the instability domains align 
themselves vertically in the (7, cos ^)— plane rather than horizontally or diagonally 
towards cos 6* = 1, 7 = 0. 

Figures Hand 121 shows the four regions in the (x, n) and (x, C) parameter planes, respectively, 
and Figs. EHSl show numerical simulations which exhibit the described behavior of the fingers. 
These figures show domains of stable and unstable behavior as a function of buoyancy x 
and rotation ( parameters for elliptic instabilities of various orders; subharmonic [n = 1), 
fundamental {n = 2), and superharmonic (n = 3,4,...). In particular, subfigures (a)-(h) 
show the results for the EB equations (To = 1) and subfigures (i)-(p) show the effects of the 
LAEB— a model for the same values of the similarity variables ( and x- 

Figure n shows the different regions for a fixed representative value of (. The horizontal 
and vertical lines which separate the different regions intersect on the parabola at the point 
(C^,C). As ( decreases, the intersection point slides down the parabola towards the origin 
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and the unstable fingers in region (I) stabilize and enter region (II). The figure, as drawn, 
indicates that the fingers corresponding to the subharmonic {n = 1) and fundamental {n = 2) 
frequencies live in region (I) for n < ^Jx^ region (IV) for ^Jx <n < ^ and in region (III) 
for n > . These fingers will be unstable in (I) and vanish in (IV) and (III). In contrast, 
the fingers for the superharmonic frequencies n > 3 live in region (II) for n < ^ region (IV) 
for <n < y/x, and finally region (I) for n > ^/x- These fingers are aligned diagonally in 
the (7, cos^)-plane for region (II) and vertically for region (IV), and then meet at the line 
7 = as parameter values enter region (I). 

We deduce two stability criteria from Figs. ^El and ()2fj|l . First, for parameter values 
satisfying [^J < y/x < [Cl ; or in original variables. 



2|n + 1| 

Tn 




2\Q+1\ 

Tn 



(29) 



where \_x\ and [x] the fioor and ceiling of x, respectively, the fiow is always stable. Here, 
fioor and ceiling refer to rounding the the nearest integer less than and greater than the 
current value, respectively. Essentially, (j^^ follows from the fact that none of the fingers 
live in region (I) for these parameter values. This criterion is analogous to the nonresonance 
stability criterion 1/2 < N/{2Q) < 2 for forced turbulence in three-dimensional rotating, 
stably stratified fiow in the Boussinesq approximation.^^ (Note that their equation for the 
buoyancy differs from ours.) For ( = \^ in Fig. HI this corresponds to 2 < y/x < 3. The 
second stability criterion is that all of the fingers will live in region (II) when 

_ _^ < + 1 < _^ for < To, (30) 

Within this window, the fiow is stable for slight perturbations in 7. This means counter 
rotation stabilizes the fiow, when the Brunt- Vaisala frequency is sufficiently low. Contrast 
this with Leblanc'ai^ stability criterion for the EB equations (Tq = 1), which says that for 
= — 1, the fiow is stable for parameter values satisfying A^^ < 1 — 7^- (In (jHUj) . 7 = 0.) 

C. Two-dimensional perturbations 



1. cos/ 



The case cos^ = 1, that is when the wave vector of the Kelvin wave is parallel to the 
axis of rotation, also can be analyzed analytically. As mentioned before, this is referred 
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to as the zero tilting vorticity because the vorticity of the elhpse-wave system is merely 
stretched rather than tilted. This follows from the fact that k = [0, 0, 1]^, which, because of 
transversality (fT^. forces a = [ai,a2,0]. The resulting equations yield b(t) = const., and 

di + (l-C + 7)a2 = 0, 

a2-(l-C-7)ai=0. (31) 

It follows that no exponential growth of the solution of ()20|) will occur for parameter values 
satisfying 

(1 - 0' - 7' > 0. (32) 
2. cos 6* = 

We perform a separate investigation for the case cos 6 = 0, that is, when the wave vector 
lies in the plane of the flow. As was the case for cos 6* = 1, the equations decouple such that 
the equations for ai 2 are independent of 03 and 6, and vice-versa. The equations for ai,2 are 
exactly those for the LANS— a model, and the equations for 03 and b can be rewritten as a 
single Schrodinger equation: 

i + ^v = 0, (33) 

where v = Ta^ and overhead dot denotes a time derivative. For the EB case, T(t) = 1, and 
the solutions are purely periodic and bounded, i.e., stable. In fact, this case is critically stable 
in the sense that any slight perturbation of the wave vector in the third dimension (cos 6 > 0) 
could result in Kelvin waves with exponentially growing amplitudes.— This corresponds to 
the fingers touching the 7— axis in the shape of a cusp, seen in subfigures (a)-(h) in Figs.OMSl 
Fig. inti-, and Fig. [7^, as predicted by Floquet theory. For the LAEB— a model, (jH^ is once 
again a Floquet problem. Comparison of numerical simulations of the full 4x4 system in 
f)20|) and that in (jH^ yield exactly the same growth rates over the same parameter ranges. 
We conclude, then, that buoyancy coupled with the model parameter Tq is the cause of this 
instability. This is illustrated in subfigures (i)-(p) of Figs. EHSl and in Fig. |H1 Numerical 
simulations suggest that these two-dimensional instabilities exist only for x < 2, even for 
extra-ordinarily large values of Tq. 
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D. Elliptic streamlines numerical simulations 

FiguresOEdisplay the instability domains for a representative set of ( and x- We reiterate 
that subfigures (a)-(h) correspond to the classic EB equations (Tq = 1), and subfigures (i)- 
(p) correspond to the LAEB— a model for parameter value Tq = 5/4. The figures are 
set side-by-side so that the reader can easily see the effects the model has on the classical 
solutions. We see the large number of instability fingers, or Arnold tongues, corresponding 
to various values of n in ()26|) . The previously described shift of fingers into and out of 
regions (I)-(IV) is also clearly seen. 

We describe the behavior of the fingers in Fig.EJ the case ( = 2 {Q = 0),m detail. For x = 
(Fig. ISt-^)) the critical angles are cos 6* = 1/2 (subharmonic) and cos 6* = 1 (fundamental), 
the latter being true for all x- The subharmonic finger lives in region (I), and the fundamental 
finger lives on the boundary between regions (I) and (II). As x increases towards unity 
(Figs. ini3;C,j,k), the subharmonic finger shifts in the parameter plane towards cos^ = 0. 
In the range 1 < x < 4 (Figs. inii,e,l,m), the finger corresponding to the subharmonic 
frequency has vanished, that is, it has entered region (IV). We conjecture that the vertically 
aligned finger in these pictures, which collapses on the line 7 = at x = 4, corresponds 
to the fundamental frequency (n = 2). In the range 4 < x < 9 (Figs. inif,g,n,o), the fingers 
corresponding to the superharmonic frequencies {n > 3) have shifted from region (II) to 
region (IV) and appear to be vertical. For x > 9 (Figs. EhiP), the fingers corresponding to 
the first superharmonic frequency (n = 3) have shifted from region (IV) to region (I) and 
meet at the prescribed critical angle (|26p along 7 = 0. Once again, the flow is unstable 
for infinitesimally small, non-zero eccentricities. Seen is the second superharmonic (n = 4) 
shifting in the plane. What cannot be seen is the fact that the finger for the subharmonic 
frequency (n = 1) has disappeared into region (III), never to return. Not shown is the fact 
that the second superharmonic frequency's fingers {n = 4) will touch on 7 = at x = 16. 
Similarly for the third superharmonic at x = 25, and so on. The fixed value of C = 2 satisfies 
()32|) for all 7, and thus the flow is always stable on the line cos 6' = 1. The vertical lines in 
Fig. ISt'-e,i-m, corresponding to the second finger, do not touch the line cos^ = 1. Rather, 
they curve sharply near that line, meeting it only at 7 = 1. This behavior is difficult to 
resolve numerically. However, we verify that this case is stable on the line cos 6* = 1 itself. 
The behavior of the fingers is similar in the presence of rotation (Figs. EJandEJ. 
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Figure m shows corresponding contour plots for the parameter value ( = 3/2. We see 
that the primary finger (subharmonic frequency) touches the 7 = line at cos 6* = 2/3 
(Figs. Note that for 1 < x < 4 (Figs. |31i,e,l,m), there are eccentric fiows for which 

the Kelvin waves are stable. This is because is satisfied for these parameter values. 
The value of ( here satisfies fl32j) only for — 1/2 < 7 < 1/2. Outside of this region, the 
fiow is unstable on the line cos^ = 1. We see that in all images, this fact holds true. Any 
visual deviation from this stability band is a consequence of the plotting program. We verify 
numerically that the band does not change width as a function of Tq or x- 

The behavior for C = 3/4, shown in Fig. is similar to the previous cases. Note that 
the main difference here is the stable band of eccentric fiows for x < 1 (Figs. int-c,i-k), 
where fl3()j) is satisfied. Furthermore, we see pairs of fingers enter region (I) as predicted-the 
primary (subharmonic) for x > 1; the secondary (fundamental) for x > 4, and the tertiary 
(first superharmonic) for x > 9 (Figs. inii-h,l-p). 

Remark. Previous work^ commented on the fact that there is a remarkable similarity 
between the effects of counter-rotation, i.e. decreasing |fi+l|, and increasing the parameter a 
on elliptic instability. This similarity is better understood by choosing ( = 2\Q+l\/{l+a'^/3^) 
as a similarity variable. Numerical simulations show that the growth rates obtained by 
maximizing over the (7, cos^)-plane are roughly constant functions of Tq for randomly 
chosen fixed values of x and ( satisfying C < 1- This is because the maximum typically occurs 
on the finger associated with the subharmonic frequency [n = 1). Again, this indicates that 
the LAEB— a model will preserve the basic characteristics of elliptic instability by properly 
adjusting the Rossby number and the Brunt- Vaisala frequency. When C > 1, however, the 
finger associated with the subharmonic frequency has vanished, and the maximum growth 
rates vary due to the presence of other frequencies. Since the model alters the shape of the 
fingers associated with these frequencies as well as the associated Lyapunov exponents for 
moderate and large eccentricities, there is no reason to expect that the growth rate remain 
constant. 
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V. ANISOTROPIC LAEB-a MODEL 

We now examine the effects of anisotropy in the third component of the velocity field v. 
We do this by altering the definition of the transported velocity v from v = (1 — a^A)u to 

v=(l-«^(C + 5,\ + e'5L))u. (34) 

Here, e measures the anisotropy of the transported velocity field < e < 1. For e = 1, 
we regain the results of the previous section. From the view of elliptic instability, the main 
difference is the change in the definition of T(t) from to 

T{t) = 1 + a^(3\kl{t) + kl{t) + e^kl{t)). (35) 

The theory for the isotropic LAEB— a model in the parameter regime I7I <C 1 remains intact 
using the similarity variables given in with 

To = l + (i-e2)cos2^^_ (36) 

Effectively, anisotropy in the third component will reduce the value of Tq (since 1 — (1 — 
e^) cos^ <1). Numerical simulations verify that there are no changes along the line 7 = 0, 
trival differences in the region I7I -C 1, and possibly noticable differences for 7^1, the 
last of which occur only for values > 1. The main distinction lies in the fact that 

each subfigure of Figs. corresponds to fixed values of f2 and A^^, whereas similar figures 
using ()23p with (jHUj) would result in different horizontal lines, that is, each value of cos^, 
corresponding to different Vt and A^ values within a subfigure. 

VI. CONCLUSIONS 

We have studied the effects of a turbulence closure model for the Euler-Boussinesq equa- 
tions on elliptic instability. We found a trade-off in the effects of turbulence, in that the 
model can increase the width of one stability band (see ((221)) while simultaneously decreas- 
ing the width of another (see ()3()|) ). This occurs because of wave number dependence is 
introduced when the smoothing operation (with lengthscale a) is applied in developing a 
mean model of turbulence. We expect that other turbulence models will show similar ten- 
dencies, although our earlier work^ shows that the magnitudes and functional forms of these 
tendencies may vary from one turbulence model to another. 
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The LAEB— a model considered here has the advantage of preserving the main features 
of classic elliptic instability analysis for circular flows, up to rescaling the stability conditions 
in terms of rotation and buoyancy frequencies by certain factors involving the nondimen- 
sional wavenumber parameter Tq, which may also contain anisotropic information as well. 
This wavenumber parameter is made nondimensional by using the lengthscale (turbulence 
correlation width) in the LAEB— a model. The similarity property and the preservation of 
structure under appropriate rescaling the frequencies by Tq has allowed a complete analysis 
of the mean effects of turbulence on the elliptic instability, when viewed as an exact non- 
linear CC solution of the LAEB— a model. This similarity property also reveals how the 
smoothing introduced in developing the turbulence model interacts with the physical effects 
of rotation and buoyancy, which are vital for understanding turbulence in geophysical flows. 
One outstanding question remains, which is something of a technical detail. Namely, we 
would still like to know the mechanism by which the turbulence model affects the resonance 
regions seen in Figs. IBE for elliptical flows of high eccenticity. These are similarity breaking 
effects of the LAEB— a turbulence model. They arise because of the introduction of the 
length scale a and cannot be removed by rescaling the other parameters in the problem. 
Anisotropy in the z— component of the velocity field amplifies these effects. These effects 
tend to occur as the shear in the flow (represented by 7) becomes more pronounced. In flows 
of such high eccenticity, coursening at the correlation length a in the turbulence model may 
be exhibiting different effects than for small eccentricity, because of the higher anisotropy of 
these flows in the plane of the strain flow. The derivation of the LAEB— a model assumes 
isotropy of the fluctuations in all coordinates, which does not hold in the case of flows with 
high eccentricity. Perhaps another version of the model-one which would allow for strong 
planar anisotropy in the statistical correlations of turbulent Lagrangian fluctuations with 
their mean trajectories- would not show such pronounced effects at high eccentricity. 

Additionally, the LAEB— a model introduces a band of unstable flows in the parameter 
regime cos 6* = 0, where the wave vector lies in the plane of the flow. This is a direct 
consequence of the presence of buoyancy. Classically, this parameter regime is critically 
stable. 

Finally, we note that Miyazaki and Adachi^i were able to extend the result of Lifschitz and 
Fabijonas^S to stratified circular flows. Namely, the stable parameter values for the Kelvin 
wave in the circular case (7 = 0) are unstable with respect to high-frequency perturbations. 
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We conjecture that the same holds true for the LAEB— a modeL 
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APPENDIX A: RECASTING EQ. (20) AS A PAIR OF SCHODINGER 

EQUATIONS 



We recast ()20|) as a coupled pair of Schodinger equations in the spirit of Refl23l for all 
base flows for which the matrix S has the form 



S 



S± O2 
Ol 



Here O2 is the zero vector in and iSj, is a 2 x 2 matrix. Incompressibility demands that 
ti{S±) = 0. 

We decompose the amplitude vector a and the wave vector k into components which are 
perpendicular and parallel to the axis of rotation: 
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where O2, and a_|_ are vectors in M^. Equations p2j) and p7 p - p9p take the form 

kT ■ + kiiaii = 0, 



-Si ■ k , 



dfh = 0, 



at(JLa^) = " b 



((T - l)Sl + S± + 2Qn 



k^kl ■ M 



a_L, 



rft(Taii) = 
dth = a\\. 



Here, 



M = {T + 1)S± + (T - i)sl + 2nn, 
^0 -1' 

1 



7^ 



Consider the transformations 

P 



Ikl 



-(kx(Ta)), 



r 



■A;||(Ta|| 
Ikl 



Ik^r ' Ik^l 

Here, (u)|| represents the third component of a vector u. The system of equations satisfied 
by these variables is 

/ 



where 





d 






dt 


1 

V7 


K = 




Q = 




Ml = 




|k± 


|k 


2 



K H Ml 
-Q -K 
-M2 -K 



(Al) 



H 



V7 

T|k|2|k^|2 



Mo 
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We eliminate the variable p from the problem, recasting the above system as a pair of second 
order Schrodinger equations for q and r (overhead dot denotes differentiation with respect 
to time): 



Q 
Q 



QH + k-K^- 



q + QMi r = 



r + Kr+[K + M1M2] r - ' ^ ' q 



Q 



+ M2H 







The results for the LANS— a model are obtained by setting A^^ = (corresponding to 
Ml = 0) and ignoring the equation for r. The results for the Navier-Stokes equations^^ are 
regained by setting T = 1. 



B. R. Fabijonas and D. D. Holm, "Mean effects of turbulence on elliptic instability," Phys. 
Rev. Lett. 90, 124501 (2003). 

^ B. R. Fabijonas and D. D. Holm, "Craik-Criminale solutions and elliptic instability in nonlinear- 
reactive closure models for turbulence," Phys. Fluids 16, 853 (2004). 

^ D. D. Holm, J. E. Marsden, and T. S. Ratiu, "Euler-Poincare models of ideal fluids with 
nonlinear dispersion," Phys. Rev. Lett. 80, 4173 (1998). 

^ D. D. Holm, "Fluctuation effects on 3d lagrangian mean and eulerian mean fluid motion," 
Physica 133D, 215 (1999). 

^ D. D. Holm, J. E. Marsden, and T. S. Ratiu, "The Euler-Poincare equations in geophysical fluid 
dynamics," in Large-Scale Atmosphere- Ocean Dynamics 2: Geometric Methods and Models, 
edited by J. Norbury and I. Roulstone, pages 251-299, Cambridge University Press, 2002. 

^ J. Leray, "Sur le mouvement d'un liquide visqueux emplissant I'espace," Acta Math. 63, 193 
(1934). 

^ Lord Kelvin, "Stability of fluid motion: rectilinear motion of viscous fluid between two parallel 

plates," Phil. Mag. 24, 188 (1887). 
^ B. J. Bayly, "Three-dimensional instability of elhptical flow," Phys. Rev. Lett. 57, 2160 (1986). 
^ A. D. D. Craik and W. O. Criminale, "Evolution of wavelike disturbances in shear flows: a class 

of exact solutions of the Navier-Stokes equations," Proc. R. Soc. London A 406, 13 (1986). 



22 



A. D. D. Craik, "A class of exact solutions in viscous incompressible magnetohydrodynamics," 
Proc. R. Soc. London A 417, 235 (1988). 

A. D. D. Craik, "The stability of unbounded two- and three-dimensional flows subject to body 
forces: some exact solutions," J. Fluid Mech. 198, 275 (1989). 

F. Waleffe, "On the three-dimensional instability of strained vortices," Phys. Fluids A 2, 76 
(1990). 

T. Miyazaki and Y. Fukumoto, "Three-dimensional instability of strained vortices in a stably 
stratified fluid," Phys. Fluids 4, 2515 (1992). 

R. R. Kerswell, "Elliptical instabilities of stratified, hydromagnetic waves," Geophys. Astrophys. 
Fluid Dynam. 71, 105 (1993). 

T. Miyazaki, "Elliptical instability in a stably stratified rotating fluid," Phys. Fluids 5, 2702 
(1993). 

S. Leblanc, "Internal wave resonances in strain flows," J. Fluid Mech. 477, 259 (2003). 
R. R. KersweU, "Elliptical instabihty," Annu. Rev. Fluid Mech. 34, 83 (2002). 
V. A. Yakubovich and V. M. Starzhinskii, Linear Differential Equations with Periodic Coeffi- 
cients, Wiley, 1967. 

C. Cambon, J. .P Benoit, L. Shao, and L. Jacquin, "Stability analysis and large eddy simulation 
of rotating turbulence with organized eddies," J. Fluid Mech. 278, 175 (1994). 
'^^ L. M. Smith and F. Waleffe, "Generation of flow large scales in forced rotating stratified 
turbulence," J. Fluid Mech. 451, 145 (2002). 

T. Miyazaki and K. Adachi, "Short-wavelength instabilities of waves in rotating stratified 
fluids," Phys. Fluids 10, 3168 (1998). 

A. Lifschitz and B. R. Fabijonas, "A new class of instabilities of rotating fluids," Phys. Fluids 
8, 2239 (1996). 

B. J. Bayly, D. D. Holm, and A. Lifschitz, "Three-dimensional stability of elliptical vortex 
columns in external strain flows," Phil. Trans. R. Soc. Lond. A 354, 1 (1996). 

24 p. N. Brown, G. D. Byrne and A. C. Hindmarsh "VODE: A Variable Coefficient ODE Solver," 
SIAM J. Sci. Stat. Comput. 10, 1038 (1989). 

25 E. Anderson et al., LAPACK Users Guide - Release 3.0, SIAM, 2000. 



23 





















(IV) 




(II) 




















(I) 
















(I),/- 






(III) 






(IV 


) 











1 4 9 16 25 



FIG. 1: The (x, n) parameter plane for the representative fixed value C = \/5- The roman numerals 
indicate different regions described in the text. Region (I) corresponds to real values of cos 6 which 
are less than unity. Parameter values in this region will yield exponential growth of a(t) for small 
values of eccentricity 7. Regions (II) and (III) correpsond to real values of cos 6 greater than unity, 
and region (IV) correpsonds to imaginary values of cos 9. Parameter values which fall into these 
regions yield windows of eccentricity values for which the disturbance is always stable. The curves 
separating different regions are the dashed curves given by the equations n = yjx (the parabola), 
X = (tti6 vertical line), and n = C, (the horizontal line). These curves separate the various types 
of elliptic instability behavior in the (x, n) parameter plane into four regions, according to Eq. H26|l 
for the critical angle. 
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FIG. 2: The (x, C) parameter plane for the representative fixed value n = 3, with the different 
regions described in the text. This figure describes the behavior of the tertiary finger for different 
values of C and x- The regions are separated by the dashed curves given by the parabola X = C^; 
the vertical line x = and the horizontal line C = n. 
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FIG. 3: Instability domain of a(t) in the classical EB equations (Figs.(a)-(h) with Tq = 1) and 
the LAEB— Q model (Figs, (i)-(p) with Tq = 5/4) for rotation parameter Q = 2 (no rotation) 
and various values of buoyancy parameter x'- (a),(i) 0; (b),(j) 1/4; (c),(k) 3/4; (d),(l) 2; (e),(m) 3; 
(f),(n) 5; (g),(o) 8; and (h),(p) 14. Regions in white represent parameter values for which the 
amplitude a is bounded in time. For all other regions, the amplitude grows exponentially in time. 
The contour lines are level lines of the Lyapunov-like growth rate with the same coloring in each 
figure. The fact that some of the fingers appear to be disjoint regions is merely an artifact of 
numerical simulation and data visualization. See the text for a description of the behavior of each 
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FIG. 4: The instability domains of a(t) in the classical EB equations (Figs.(a)-(h) with Tq = 1) 
and the LAEB— a model (Figs.(i)-(p) with Tq = 5/4) for rotation parameter ( = 3/2 and the same 
values of buoyancy parameter x as in Fig. |21 (a),(i) 0; (b),(j) 1/4; (c),(k) 3/4; (d),(l) 2; (e),(m) 3; 
(f),(n) 5; (g),(o) 8; and (h),(p) 14. Note the presence of a stable band in (d)-(e) and (l)-(m), 
corresponding to 1 < ^/x < 2 (see (EHI))- 
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FIG. 5: The instability domains of a(t) in the classical EB equations (Figs.(a)-(h) with Tq = 1) 
and the LAEB— a model (Figs.(i)-(p) with Tq = 5/4) for rotation parameter = 3/4 and the same 
values of buoyancy parameter x as in Fig. |21 (a),(i) 0; (b),(j) 1/4; (c),(k) 3/4; (d),(l) 2; (e),(m) 3; 
(f),(n) 5; (g),(o) 8; and (h),(p) 14. Note that since C < 1, we have a stable band of flows in (a)-(c) 
and (i)-(k), where x < 1 (see H3U|) ). 
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FIG. 6: Instability domain of a{t) for rotation parameter C = 3, buoyancy parameter x = 1/4, and 
various model parameter Tq: (a) 1 (the EB equations), (b) 3/2, (c) 2, (d) 11/4, (e) 4, and (f) 6. 
The behavior of the instability regions is essentially the same in all figures for I7I <C 1. 
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FIG. 7: Instability domain of a{t) for rotation parameter ^ = 3/2, buoyancy parameter x = 5, and 
various model parameter Tq: (a) 1 (the EB equations), (b) 3/2, (c) 2, (d) 11/4, (e) 4, and (f) 6. 
Again, the behavior of the instability regions is essentially the same in all figures for I7I <C 1. 
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FIG. 8: The growth rate for the case cos6 = for x = 3/4 and ^ = 3, and for values of Tq in 
the range 1 < Tq < 4. For Tq = 1, the Lyapunov exponent is identically zero. The exponent 
is immediately non-zero for Tq > 1, and shifts to the left as the parameter increases. For larger 
values of Tq, the growth rates approach an asymptotic limit, which is slightly to the left of the 
curve corresponding to Tq = 4 shown in the figure. 
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FIG. 9: The growth rate maximzed over the {cos 9,^) plane as a function of x for various ^ and 
Tq: (a-b) To = 1, (c-d) Tq = 3/2, (e-f) Tq = 4. The figures in the left column correspond to 

values satisfying < C ^ 1; increasing from bottom to top. Similarly, the figures in the right 
column correspond to (" values satisfying <^ > 1, increasing from top to bottom. Notice that as 

— > oo, the growth rates approach an asymptotic limit for all Tq. This is in accordance with the 
Taylor-Proudman theorem. Notice further the development of a "bump" for Tq > 1, when C = 
(the bottom-most curve in the figures in the left column) and as ^ oo (the bottom-most curves 
in the figures in the right column). 



